COMPRESSIVE IMAGING OF SUBWAVELENGTH STRUCTURES 



ALBERT C. FANNJIANG 

Abstract. The problem of imaging extended targets (sources or scatterers) is formulated in the 
framework of compressed sensing with emphasis on subwavelength resolution. 

The proposed formulation of the problems of inverse source/scattering is essentially exact and 
leads to the random partial Fourier measurement matrix. In the case of square-integrable targets, 
the proposed sampling scheme in the Littlewood-Paley wavelet basis block-diagonalizes the scat- 
tering matrix with each block in the form of random partial Fourier matrix corresponding to each 
dyadic scale of the target. 

The resolution issue is analyzed from two perspectives: stability and the signal-to-noise ratio 
(SNR). The subwavelength modes are shown to be typically unstable. The stability in the subwave- 
length modes requires additional techniques such as near-field measurement or illumination. The 
number of the stable modes typically increases as the negative d-th (the dimension of the target) 
power of the distance between the target and the sensors/source. The resolution limit is shown to 
be inversely proportional to the SNR in the high SNR limit. 

Numerical simulations are provided to validate the theoretical predictions. 



1. Introduction 

The Compressed Sensing paradigm is supplying a fresh perspective on the imaging problems, 
including source inversion and inverse scattering (see the extensive literature in [8]) and we aim to 
analyze these classical problems from the perspective of compressed sensing. 

The main purposes of the paper are two folds: (i) to formulate the problems of inverse source 
and scattering for extended targets in the framework of compressed sensing and analyze it by 
the fundamental results on the random Fourier measurements [H El El [3 EH |22] and (ii) to use 
the compressed sensing solution thus obtained as a starting point to investigate the problem of 
subwavelength resolution in the presence of noise, in particular, to discuss the idea of extracting 
subwavelength information by near-field measurement and illumination, common in nano-optics 

m- 

The important feature of our formulation of the imaging problems is that it is essentially exact 
(without the paraxial approximation as used in [TH]) and is for extended (periodic or non-periodic) 
targets. As such, our formulation provides an instructive example for assessing the power of the 
compressed sensing techniques when applied to the physical problems of imaging. In the case 
of inverse scattering for square-integrable targets, we use the Littlewood-Paley wavelet basis and 
propose a sampling scheme by which the scattering matrix can be block-diagonalized with each 
block in the form of random partial Fourier matrix corresponding to each dyadic scale present in 
the target structure. By this approach we can image the extended target scale-by-scale. 

We analyze the resolution issue from two perspectives. First we define the stably recoverable (or 
stable for short) modes in terms of the noise amplification factor. We show that the subwavelength 
modes are typically unstable. To achieve stability in the subwavelength reconstruction we need 
additional techniques such as near-field measurement or illumination with which the number of the 



The research is partially supported by the NSF grant DMS - 0908535. 

1 



stable modes typically increases as the negative d-th (the dimension of the target) power of the 
distance between the target and the sensors/source. Second, we take into account the signal-to- noise 
ratio (SNR) in the resolution analysis and show that the resolution limit is inversely proportional 
to the SNR in the high SNR Umit. 

The paper is organized as follows. We treat the case of periodic targets using the Fourier basis, 
first for the inverse scattering in Section |2] and then for the inverse (Born) scattering problem in 
Section 3 , We present numerical results confirming the theoretical predictions in Section |4} In 
we analyze the case of square-integrable targets using the Littlewood-Paley basis and 



Section 



point out the extension to three dimensions. We conclude in Section [6} 

2. Source inversion 

First we consider imaging of an extended source (the target) in two dimensional (z, x)-plane. We 
assume that the target is located in the line z = zq. 

The target is represented by the variable source amplitude cr(x) which is assumed to be periodic 
with period L and admits the Fourier expansion 



(1) a{x) = ^ cjfce" 



k=—oo 

We consider the class of band-limited functions with = except for 

(2) k £[-{m-l)/2,{m-l)/2] 

in ([T| where m is an odd integer so that there are at most m relevant Fourier modes. 
The wave propagation in the free space is governed by the Helmholtz equation 

Au + up'u = 0. 

where oo = 2tt/\ is the frequency assuming the wave speed is unity. In two dimensions, the Green 
function G{y) is 

G{v)='-^H^^\u\v\), r = {z,x) 

where Hq^^ is the zeroth order Hankel function of the first kind. G can be expressed by the 
Sommerfeld integral representation 

(3) G'(r) = J_ f ^M\zma)+xa) da 



4vr J (3{a) 
where 



(4) /3(a) 



iV — 1, IqI > 1 



[T]. The integrand in ([s]) with real- valued (3 (i.e. \a\ < 1) corresponds to the homogeneous wave 
and that with imaginary-valued /3 (i.e. \a\ > 1) corresponds to the evanescent (inhomogeneous) 
wave which has an exponential-decay factor g-'^klVa-i_ 

The signal arriving at the sensor located at (0,x) is given by 

(5) [ G{zo,x-x')a{x')dx' = ^ ^ ^g^'^^o^^e^'^"^^' 

J 2uj ^ (3k 

by ^ where 

k\ 

(6) "'^ ^ IT' Pk = P{ak)- 
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By the subwavelength structure of the target is encoded in with |A;|A > L which are conveyed 
by the evanescent waves. 

Let (0, Xj),Xj = (,jL, j = 1, n be the coordinates of the n sensors in the hne z = 0. To set the 
problem in the framework of compressed sensing we define the signals received by the n sensors as 
the measurement vector Y and set the target vector X = (X^) 



To avoid a vanishing denominator in Q, we must set L/A N so that Pk / 0. This leads to the 
form 

(8) Y = ^X 
where the sensing matrix $ = [^jk] G C"^™ has the entries 

(9) '^jk = ^e^-^^^ = ^e^2.fcc„ j- = i,...,^, k = l,...,m. 



Hence the imaging procedure is split into two stages: the first is to "invert" ^ and solve for X, and 
the second is to deduce o from X. The sensing matrix ([9| is designed to mimic the random partial 
Fourier matrix which has the restricted isometry property and is amenable to the compressed 
sensing techniques (see below). Note that the formulation is essentially exact and no far- field 
approximation is made. 

The main thrust of compressive sensing is that under suitable conditions the inversion can be 
achieved as the ^-'^-minimization 

min||X||i, subject to Y = ^X 

which also goes by the name of Basis Pursuit (BP) and can be solved by linear programming. BP 
was first introduced empirically in seismology by Claerbout and Muir and later studied mathemat- 
ically by Donoho and others |2l |9l [13] . 

A fundamental notion in compressed sensing under which BP yields the unique exact solution 
is the restrictive isometry property due to Candes and Tao |6]. Precisely, let the sparsity s of 
the target vector be the number of nonzero components of X and define the restricted isometry 
constant 5s to be the smallest positive number such that the inequality 

{l-6s)\\Zg<\\^Zg<{l + 5s)\\Z\\l 

holds for all Z G of sparsity at most s. 

For the target vector X let Xs denote the best s-sparse approximation of X in the sense of 
L^-norm, i.e. 

X'^ = argmin \\Z — X\\i, s.t. ||.^||o < s 

where \\Z\\q denotes the number of nonzero components, called the sparsity, of Z. Clearly, X^ 
consists of the s largest components of X. 

Now we state the fundamental result of the RIP approach [3j which is an improvement of the 
results of HE]. 

Theorem 1. [3 Suppose the restricted isometry constant of ^ satisfies the inequality 



(10) 



(52« < \/2 - 1 
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Then the solution by BP satisfies 

(11) ll^*-^lli < Co\\X-X'\\i 

(12) \\X^-X\\2 < Cos"^/^||X-X^||i 

for some constant Cq. In particular, if X is s-sparse, then the recovery is exact. 

Remark 1. Greedy algorithms have significantly lower computational complexity than linear pro- 
gramming and have provable performance under various conditions. For example under the condi- 
tion < 0.06 the Subspace Pursuit (SP) algorithm is guaranteed to exactly recover X via a finite 
number of iterations [14j. See [19j for a closely related algorithm (CoSaMP). 

Theorem 2. [21] Let = 1, ...,n be independent, uniform random variables in [0, 1]. Suppose 

(13) > Ct^^Vln^rlnmln-, ?? G (0, 1) 

In n rj 

for a given sparsity r where C is an absolute constant. Then the restricted isometry constant of the 
matrix for random Fourier measurements satisfies 

6r<d 

with probability at least 1 — rj. 

See [H [3 [22] for the case when belong to the discrete subset of [0, 1] of equal spacing 1/ y/m. 

As a consequence of Theorem [T] and |2] the target vector X can be determined by BP and then 
the source amplitude cr(x) can be reconstructed exactly from ([7]), including all the subwavelength 
structures. In other words, in the absence of noise there is essentially no limitation to the resolving 
power of the compressed sensing technique, subwavelength or not, as long as sufficient number of 
measurements are made. 

When noise is present, however, the performance of the above approach may be severely limited, 
especially in the recoverability of subwavelength information. Consider the standard model of 
additive noise 

(14) Y'' = ^X + E 

where ||£^||2 = £ > is the size of the noise and the associated relaxation scheme 

(15) min||X||i, subject to - ^X\\2 < e. 

The next result is a restatement of the result of [3] after applying Theorem[2]with r = 2s, 5 < \f2 — \. 



Theorem 3. Let X^ be the solution to {'15). Then under the assumptions of Theorem ^ we have 

(16) \\X^ - Xy < Cqs'^^^WX - X^lli + Cie 
with probability at least 1 — rj where Cq and Ci are constants. 

See also |12| l23 j for stability result under the condition of incoherence. 

Inverting the relationship ([T]) with small error in the target vector X produces a mildly amplified 
error for those ak such that 

(17) |gi^20A| > g-27r 

but significantly amplified error otherwise. Here the transition is not clear-cut, however. The choice 
of the noise amplification threshold e^'^ as the stability criterion is convenient but arbitrary; any 
constant less than one will serve our purpose. 
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The stable (or stably recoverable) modes are those corresponding to |afc| < 1 as well as {a^l > 1 
such that 



(18) 

or equivalently 
(19) 



k < L\ \-^ + z, 







Hence reducing the distance zq between the sensor array and the target can effectively enlarge the 
number of stable modes and more of the subwavelength modes become stably recoverable as zq 
decreases below wavelength. This is the idea behind the near-field imaging systems such as the 
scanning microscopy. 

The above reconstruction essentially consists of two stages: Stage 1 involves the compressed 
sensing techniques which is always stable (Theorem |3]) and Stage 2 is a simple inversion of the 
diagonal matrix 



(20) 



diag 



cf. ([T]), which is stable if and only if (18) holds for all k such that (Tk 7^ 0. When (18) is indeed 
violated, brute force inversion of the matrix (20) would lead to an enormous error and hence a 
regularization is called for. In such a case, we can use the well-known Tikhonov regularization to 
invert (|20). 



The analysis from (17) to (19) focuses on the stability issue based on the noise amplification 
factor. On the other hand, it is well accepted that the resolution of any imaging system should 
be a function of the signal-to-noise-ratio (see flT] and references therein) . Let us now analyze the 
resolution from this perspective. 



(21) 



Let X'' = {Xf,) be the solution of the convex relaxation scheme (15) and let 

iJn 



For simplicity, let us assume that the sparsity of X is at most s and therefore by Theorem [3] we 
have 

(22) \\X^-X\\2<Cie 

where the constant Ci is well-behaved as demonstrated in the numerical results of Section |4j 
We say that the Fourier mode is resolved if the following inequality holds: 

1 



1 



< 



K 



where the constant K is sufficiently larger than unity. Using (22) and (21) we obtain 



CFl. 



1 



< 



n 



Xh 



(23) 

The expression 



< 2Ciuj\l3k\ 



e 



SNRfe 



Ju)ZoPk 



\akWn/e 



has the meaning of the signal-to-noise ratio for the k-th Fourier mode and we shall call it as such. 
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Assuming that the right hand side of (23) is at most 1/K we obtain from (23) the inequahty 

1*1 < 

and consequently 

/ . X 9\ -1/2 

L I f SNRfc ^ 



The resolution limit is then obtained by minimizing the left hand side of (25). For simplicity, 
considering the case that SNR^ = SNR is independent of k for all nonzero modes dk we obtain the 
resolution limit 

/ . ^ 9\ -1/2 

^ 1 / SNR ^ 



which essentially says that the resolution is inversely proportional to SNR in the high SNR limit. 

Next let us discuss a different mechanism of superresolution available in the context of inverse 
scattering. 

3. Inverse scattering 

Consider a periodic scatterer with scattering amplitude a admitting the representation ([T|. For 
simplicity, we use the Born scattering model [1] under which the scattered field at z = is given 
by 

(27) u'^(0,x)= / G{zq^x — x')a{x')u^{zQ,x')dx 



where u^{z,x) is the incident field. Now with the normally incident plane wave v}{z,x) = e**^^, eq. 
(27) is essentially reduced to 

Consider the obliquely incident plane wave u'(z,x) = where /5 is related to a as 

in (|6| and zi, assumed larger than zq. Set 

In the standard setting the illumination field is a homogeneous wave with \a\ < 1. The evanescent 
illumination |a| > 1 will also be considered here. In such case, zi is the z-coordinate of the 
illumination source. 

The same calculation as before now leads to 

where instead of ^ we have 



{k + q)X I J^-al \ak\ < 1 

«fc = 7 ' Pk= < I— 

^ I iJal - 1, lofcl > 1 



Define the target vector X^") = [x]^^) G C"" as 



k 



(28) X]^> = — a. 



~2^k 



and proceed as before. Note that e*'^^^^ ^"^^ is a constant factor determined solely by the illumi- 
nation source (distance and angle). To avoid a vanishing denominator we require 

(A; + g)A 



(29) 



L 



/ 1, yk. 



Theorem [l|and|2| are applicable to the shifted Fourier matrix = [^^^^] with 



(<?) 
jk 



^i2-K{k+q)ij 



n 



where are independent, uniform random variables in [0, 1]. 

The same stability analysis as in Section [2] implies that the stable modes in the Born scattering 
satisfy the constraint 



(30) 

which implies that for each q 
(31) -q 



\k + q\X 
L 



< 



'L2 L2 , 



'1 + 



+ 



A2 



'L2 L2 



A2 



+ 



Let us maximize the range of stable modes (31) under the constraint |e 



iuj{zi-zo)P\ 



> e 



-2n 



by 



considering two different illumination sources: homogeneous and evanescent wave sources. If the 
target is illuminated by incident homogeneous waves with l^l < L/X, then the stable modes are 



(32) 




The same resolution analysis as before leads to 



(33) 



< 



2w|/3fc| 



n 



-iLu{zi-zo)l3 



The appropriate definition of the signal-to-noise ratio in this case is 



SNRi 



n\ak\ 



JuizofSk 



Juj{zi-zo)l3 



after taking into account the intensity of illumination. Setting the right hand side of (33) to be at 
most ^jK then leads to (24) with given by (30). which in turn yields the resolution limit 

-1 



(34) 



+ 



1 



+ 



SNR 



A2 X^TlCxK 



On the other hand, if the incident wave is evanescent and subject to the constraint |e*'^'^^i ^"^^1 > 
€"2^^^ which implies 



(35) 




then the stable modes according to (31) satisfy 

, , ,,, ,^ ^ 

(36) 




1=1001, n-35, s-11, z,=1.0, lambda=0.104720 




m=1001, n-35, s-11, z =0.1. lambda=0. 104720 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 1. In source inversion with far-field measurement (left, zq = l,n = 35), 
the two subwavelength modes k = ±11 cause significant errors. With near- field 
measurement (right, zq = 0.1, n = 35) the reconstruction is nearly perfect for target 
with 20 subwavelength modes. The thin-red curve is the original profile and the 
thick-green curve is the reconstructed profile. The two curves coincide in the plot 
on the right. 



To maximize the range of stable modes limited by (36) we choose > 0, which is sufficiently 



close to the right hand side of (35) and satisfies (29), and illuminate the target by a few evanescent 
waves with q £ [—Q*,Q*]- The union of their respective stable modes (31) is the totality of stable 
modes. In optics, the evanescent illumination can be produced physically by, for example, the total 
internal reflection [20j. 

Using the same evanescent illumination procedure we can achieve the resolution limit 

-1 

/ / 1 1 / 1 / CMt> \ 2 \ 

(37) 

which indicates two ways of improving resolution without increasing the probe frequency: reducing 
the distance between the plane-wave source and the target, and increasing the signal-to-noise ratio. 




4. Numerical results 

Numerical methods for convex optimization (such as the Matlab program cvx designed by M. 
Grant and S. Boyd) are not exact; they compute their results to within a predefined numerical 
precision or tolerance which is acceptable for most applications. However, for our calculations 
involving subwavelength modes, the numerical error tends to be amplified enormously and hence 
severely spoils the reconstruction. 

On the other hand, because of the simplicity of the greedy algorithms, the numerical error can 
be much reduced and the quality of recovery much improved. Of course, the sparsity constraints 
of the greedy algorithms tend to be more severe than that of the convex relaxation. Since the 
sparsity constraint is not the main focus of our study, we use the Subspace Pursuit in the numerical 
calculations. 

In our numerical simulations, L = 1, A = 7r/30, m = 1001 so the subwavelength mode cutoff is 
at about \k\ = 10. The large m is selected to demonstrate the compressed sensing step. 

The numerical results for source inversion are shown in Figure [T] First we image a periodic 
source with 11 modes, including one subwavelength mode A; = 11 using far-field measurement 
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Figure 3. Two independent runs (left with n = 70 right with n = 82) with far- field 
measurement (zq = 1) and near-field illumination {zi — zq = 0.01). The red curve is 
the original profile and the blue curve is the reconstructed profile. Except for a few 
spots, they coincide with each other. 



{zq = l,n = 35). The main source of noise in this simulation is the roundoff error. As evident in 
the left plot of Figure [T] the recovery is not accurate. Indeed, the subwavelength mode is entirely 
missing in the reconstruction. Reducing the distance between the target and the sensors {zq = 0.1) 
enables accurate reconstruction (right plot in Figure [T]). When additional random noise is added, 
the quality of recovery with near-field measurement deteriorates in proportion to the amount of 
noise (Figure |2]). With about 5% random noise, the imaging result with near-field measurement 
(Figure [2] right) is comparable to that of far- field measurement without added noise (Figure [T| 
left). 

Figure [3] shows the results for inverse Born scattering by the far-field measurement {zq = 1) and 
the near-field illumination (zi — zq = 0.01) without additional noise. In this case it suffices to run 
the procedure for two incident modes with q = ±90. 



5. Extension and generalizations 
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5.1. Non-periodic targets. Let us consider a non-periodic target represented by a square-integrable 
function a. 

To represent a we use the Littlewood-Paley basis 

(27r)"V2^ vr < 1^1 < 27r 
0, otherwise 



(38) m 



or 

(39) ip{x) = (7rx)~"^(sin (27rx) — sin (ttx)). 
Then the following set of functions 

(40) i^p,,{x)=2-P/^i^{2-Px-q), p,qGZ 

forms an orthonormal wavelet basis in L^(M) [11 . Expanding the target profile ct(x) in the 
Littlewood-Paley basis {ipp,q} we have 

(41) a{x) = (Tp,qijp,q{x). 

p,q&Z 

The scattered field has the far-field asymptotic 

jw|r| 

(42) u%r) = j-^{A{i,d,u;)+0{\r\-^)), f- = r/|r| 

where the scattering amplitude A (at frequency lo and the sampling direction r with the incident 
field = e**^'"") can be calculated according to 

(43) A{^,d,uj) = ^ /'drV(r')e^'"'^■'■'e-^'^'■'•^ d = (a,/?), + /?2 = i 

47r 7 



In inverse scattering theory, the scattering amplitudes are measured in the directions ffc = 
(afc,/3fc) = {cos6k,sm9k), k = l,...n and denoted by y = (Yfc) G C". By (43) and (41) we have 

(44) yk = -^ V 2f/%,,e^-("-"'=)2^VMa-afe)2^'), k = l,...,n. 
Let 

p-i 

I = 'Y{2mj + l) + q, \q\<mp, \p\<p*, 
j=-p* 
p'-i 

k = 'Y {2nj + 1) + q , \q'\<np, \p'\<p* 

for some mp,np,p^: € N be the column and row indices of the sensing matrix, respectively. Set 
m = Eh<p.(2"1p + 1) and n = Eh<p.(2'^p + !)• 
Suppose that 

(45) a;2~P* > 27r 

and consider, for simplicity, the normal incident field a = 0. Let Cp',g' be independent, uniform 
random variables on [—1, 1] and let 

(46) ak- — -^ ^ ^ ^P''^' ^ 



uj2p' 1-1 + Cp',g', Cp',g' G [-1,0] 
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which hes in [—1, 1] by the assumption (45). Let the sensing matrix elements be 

(47) = ^=i==^Mfc2f)e— "'=2^". 

^y2np + 1 

We claim that ^k^i = for p ^ p'. This is evident from the following calculation 

^ ' \-l + Cp',g', Cp'.g' e [-1,0] 

whose right hand side is outside the support of tjj for p ^ p' ■ 
For every p = p' , 

(49) ^k,l = -j=L=e~'<v',."i, y\ < Up, \q\ < rup 

^2np + 1 



which constitute the random partial Fourier matrix. In other words, under the assumption (45) 
the sensing matrix $ = [$fc,/] G ^^x*^ block-diagonal with each block (indexed by p) in the form 
(49]). 

Let X = [X]) be the target vector with 

2 p— 1 

^/ = ^^^(2np + l)2Pap,g, / = + 1) + 

i=-p. 

We can then express the measurement vector Y = (Y^.) E as in (|8|. 

Each block represented by the submatrix (49) can be treated as in the periodic case. For each 
Vi \v\ ^ P*) let Sp be the sparsity of the target vector for the scale 2^ and suppose that the frequency 

used to probe the scale 2"^ satisfies 



(50) u)p2^ > 2-K. 
Under these assumptions, Theorem [2] yields the sufficient condition 

(51) -4?^^^ > C^-^Spln^Splnmln-, r/ G (0, 1) 
for the RIP 

to hold with probability at least 1 — rj. 

The assumption (50) means that the wavelength is at least as small as the scale being probed. 
Therefore, this imaging method does not possess the subwavelength resolving power. 

Figure |4] shows the result of reconstruction with the Littlewood-Paley basis and the following pa- 
rameters: p^ = 2,mp = 100, Vp; for p = -2,-1,0,1,2, Sp = 12, 24, 13, 24, 23, = 36,64,36,64,64, 
e = 0.5134, 1.4849, 0.6520, 1.0274, 1.3681 equivalent to the percentages of noise = 0.0500, 0.0649, 0.0494, 0.0640, 0.06 
The resulting reconstruction errors divided by the noises are 0.1941,0.1240,0.2498,0.2361,0.1502, 
demonstrating the stability of the recovery. 

5.2. Three dimensions. The present framework can be easily extended to three dimensions by 
using the plane-wave representation 

(52) -— = — / exp[iu;{ax + I3y + j\z\)\, r = {x,y,z) 

\r\ 271 J 7 
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Figure 4. Imaging of a non-periodic target in the presence of noise: The red-sohd 
curve (top) is the exact profile and the blue * shows the reconstructed profile; the 
bottom plot shows the errors in reconstruction in the presence of noise. 

where 

7 = + /32 - 1, a^ + p"^ >l 

for the Green function [l]. 

6. Conclusion 

In this note, we have analyzed the problem of imaging extended targets in the perspective of 
Compressed Sensing which provides assurance of stable reconstruction of a target composed of 
sparse Fourier modes with the similar number of measurements, modulo a poly-logarithmic factor, 
in the presence of noise. 

We have shown that the number of stably recoverable modes grows as the negative d-th power 
of the distance between the target and the sensors /source. Hence the stability of reconstructing 
subwavelength modes requires the distance to be less than the wavelength. On the other hand, we 
have also shown that the resolution limit is inversely proportional to the SNR in the high SNR 
limit. As a consequence, the subwavelength modes can be recovered at sufficiently high SNR as 
well as by placing the sensors at a subwavelength distance from the target. 

However, it remains to be seen if these results about stability and resolution can be extended 
to the case of square-integrable targets even though our sampling scheme block-diagonalizes the 
corresponding scattering matrix according to each dyadic scale. 
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Finally we note that the compressive imaging theory in the remote sensing regime for discrete 
point targets has been recently developed in [151 HI ITS] - 

and Wenjing Liao (Figm'e 
|4j) for preparing the figures. I am grateful to the anonymous referee for the helpful comments for 
improving the manuscript. 
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